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Abstract 

This  research  investigates  the  generation,  display,  and  interpretation  of  three-dimensional 
(3-D)  Synthetic  Aperture  Radar  images.  Three-dimensional  assumes  that  the  data  collected 
consists  of  one  temporal  dimension  and  two  orthogonal  angular  dimensions.  From  this  data,  a 
three  dimensional  reflectivity  map,  or  3-D  image,  of  a  target  can  be  constructed.  This  reflec¬ 
tivity  map  can  then  be  compared  to  the  known  structure  of  a  target  to  gain  more  insight  into 
the  scattering  mechanisms  of  the  target.  Until  now  research  in  this  area  has  been  limited  to 
two-dimensional  imaging.  With  the  added  third  dimension,  a  more  detailed  target  scattering 
map  is  possible. 

First,  a  detailed  review  of  radar  imaging  theory  is  conducted.  This  review  starts 
with  basic  one-dimensional  down-range  profiles  and  Synthetic  Aperture  Radar  theory.  Two- 
dimensional  imaging  theory  is  then  presented  and  developed  through  two-dimensional  Filtered 
Back-Projection.  This  two-dimensional  theory  discussion  is  used  as  a  stepping  stone  for  the 
derivation  of  the  three-dimensional  Filtered  Back-Projection  algorithm.  Finally,  this  algorithm 
is  converted  to  a  discretized  algorithm  for  use  on  a  computer. 

This  thesis  effort  applies  this  three-dimensional  algorithm  on  actual  radar  data  measured 
on  an  I  scale  model  of  a  C-29  aircraft.  In  the  process  of  implementing  this  algorithm  a  few 
problems  were  encountered.  Mainly,  due  to  the  vast  quantity  of  data  to  be  processed,  the 
amount  of  memory  available  in  terms  of  hard  drive  space  as  well  as  RAM  became  a  limitation. 
Problems  arose  from  the  need  to  upsample  the  projection  data  in  order  to  get  sufficient  image 
quality.  Working  within  these  memory  constraints,  three-dimensional  images  were  produced. 

The  results  demonstrate  the  ability  to  produce  three-dimensional  images  given  three- 
dimensional  radar  data.  Two-dimensional  slices  of  the  three-dimensional  image  as  well  as 
three-dimensional  isosurfaces  are  compared  to  the  physical  properties  of  the  target. 
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/.  Introduction 

1.1  Problem 

The  United  States  Air  Force  is  actively  trying  to  develop  technology  to  produce  high 
resolution  spatial  maps  from  radar  data  taken  on  various  targets.  These  maps  can  then  be 
analyzed  for  target  identification  purposes.  A  couple  of  recent  incidents  reinforce  the  USAF’s 
interest  in  this  technology.  Several  years  ago,  the  United  States  Navy  accidently  shot  down 
an  Iranian  747,  mistaking  it  for  a  hostile  fighter  aircraft.  More  recently,  the  USAF  shot 
down  two  United  States  Army  Blackhawk  helicopters  mistaking  them  for  much  larger  Iraqi 
Hind  helicopters.  The  obvious  structural  differences  between  these  mistaken  identities  and 
the  actual  targets  could  have  been  discerned  had  this  radar  technology  been  operationally 
available. 

To  this  point,  SAR  imaging  research  has  been  limited  to  two-dimensions.  With  the 
added  third-dimension,  more  information  can  be  obtained  about  a  target,  ultimately  allowing 
a  more  accurate  target  identification. 

This  research  investigates  the  generation,  display,  and  interpretation  of  three-dimension 
(3-D)  Synthetic  Aperture  Radar  images.  Three-dimension  assumes  that  the  data  collected 
consists  of  one  temporal  dimension  and  two  orthogonal  angular  dimensions.  From  this  data, 
a  three  dimensional  reflectivity  map,  or  3-D  image,  of  a  target  can  be  constructed. 

Advanced  radar  ranges  can  collect  and  process  data  of  this  type  efficiently.  This  research 
develops  an  efficient  imaging  algorithm  to  process  this  data.  Also,  since  the  resulting  image 
is  a  function  of  three  dimensions,  methods  to  properly  display  the  images,  or  slices  of  the 
images,  will  also  be  developed. 
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1 .2  Radar  Imaging  Background 

A  thorough  knowledge  of  the  current  two-dimensional  (2-D)  research  is  required  in 
order  to  the  develop  an  extension  to  3-D.  This  summary  provides  an  overview  of  past  research 
done  in  this  area. 

1 .2.1  Parametric  Reconstruction  Methods. 

In  many  applications  of  computerized  imaging,  such  as  in  electron  microscopy,  astron¬ 
omy,  geophysical  exploration,  non-destructive  evaluation,  and  others,  it  is  not  possible  to 
collect  a  complete  set  of  angular-projections.  Standard  Fourier  based  reconstruction  methods 
perform  poorly  in  the  case  of  limited  projection  data  due  to  poor  resolution.  Parametric  or 
iterative  algorithms,  on  the  other  hand,  are  capable  of  producing  high-quality  images  with 
limited  data  sets  since  they  automatically  take  into  account  the  missing  data,  as  well  as  use 
a-priori  information  about  the  target  (3). 

The  maximum  likelihood  method  (MLM)  of  spectral  estimation  is  one  parametric 
method  applied  to  the  tomographic  problem.  This  method  has  been  found  to  provide  high 
resolution,  but  does  not  exist  for  all  data  sets  (11).  Two-dimensional  auto-regressive  (AR) 
modeling  algorithms  also  exist  but  have  proved  to  be  computationally  expensive  (6).  The 
Prony  method  coupled  with  total  least  squares  (TLS)  has  been  used  successfully  to  estimate 
frequencies  from  noisy  data  in  one-dimension.  This  method  was  also  extended  to  2-D  (11). 

Delaney  and  Bressler  (3)  develop  a  series-expansion  approach  that  is  used  in  a  fast 
and  iterative  tomographic  reconstruction  algorithm.  This  algorithm  is  applicable  for  parallel- 
ray  projections  that  have  been  collected  at  a  finite  number  of  view  angles.  A  conjugate 
gradient  algorithm  is  used  to  minimize  a  regularized  least-squares  criterion.  The  main  step  in 
each  iteration  is  equivalent  to  a  2-D  discrete  convolution,  which  can  be  cheaply  and  exactly 
implemented  via  the  Fast  Fourier  Transform  (FFT)  (3). 

These  iterative  or  parametric  techniques  have  been  shown  to  have  superior  resolution 
capabilities  but  have  also  shown  to  be  computationally  expensive  (2).  This  computational 


2 


cost  is  a  major  drawback  since  just  one  iteration  typically  requires  more  computation  than  an 
entire  typical  Fourier  based  method  (3). 

1 .2 .2  Fourier  Reconstruction  Methods. 

In  Fourier  based  methods,  the  radar  returns  from  a  target  can  be  considered  as  a  sample 
of  the  2-D  Fourier  spectrum  of  the  target’s  reflectivity.  The  imaging  problem  becomes  one  of 
reconstructing  the  spatial  reflectivity  of  the  target  given  these  discrete  samples  of  its  Fourier 
spectrum  (1). 

One  problem  that  arises  in  the  use  of  Fourier  based  methods  is  that  radar  data  is  not 
usually  collected  in  a  uniform  rectangular  grid  in  the  spectral  domain  since  the  rotation  of  the 
radar  around  the  target  introduces  polar  aspects.  This  property  causes  distortion  if  the  standard 
2-D  Inverse  Fast  Fourier  Transform  (IFFT)  is  used  to  approximate  the  Fourier  Transform, 
especially  when  data  is  obtained  over  large  look  angles.  The  2-D  IFFT  is  considered  to  yield 
the  poorest  results  in  terms  of  image  quality,  but  is  extremely  fast  computationally  ( 1 ) .  In  order 
to  eliminate  the  distortion  a  focusing  procedure  must  be  used  (8).  A  few  focusing  methods 
are  now  described. 

1.2. 2.1  Direct  Focused  Imaging. 

Focused  imaging  utilizes  a  process  which  applies  exact  phase  corrections  to  the  sampled 
reflectivity  data.  Two  direct  focused  methods  are  studied  in  (1);  the  first  is  simply  an  inverse 
Fourier  transform  from  polar  coordinates  to  rectangular  coordinates.  The  second,  originally 
derived  by  Mensa  (8),  is  a  similar  method,  but  the  phase  correcting  data  points  are  not  properly 
weighted  in  the  final  sum  (1).  This  method  still  gives  very  accurate  results,  especially  as  the 
ratio  of  center  frequency  to  bandwidth  increases.  The  fact  that  an  FFT  algorithm  cannot  be 
used  causes  these  direct  focusing  methods  to  be  computationally  expensive  but  provides  the 
most  accurate  image  possible. 


3 


Technique 

Focused  (both  types) 

552 

Unfocused 

1.74 

Focused  via  Interpolation 

2044 

Filtered  Back-projection 

217 

Table  1.  Computational  Speed  of  Imaging  Techniques. 

1.2. 2. 2  Focused  Imaging  via  Interpolation. 

In  order  to  use  the  FFT  to  increase  computational  speed,  the  polar  data  is  mapped  onto 
a  rectangular  grid.  This  mapping  is  done  by  up-sampling  the  data  set  and  then  picking  the 
closest  data  point  to  the  required  rectangular  coordinates.  This  process  obviously  produces 
some  inaccuracies  and  has  been  shown  to  be  very  computationally  expensive  (1). 

1.2. 2. 3  Filtered  Back-Projection. 

An  estimated  image  can  be  constructed  by  summing,  or  superimposing,  1-D  projections 
from  all  available  angles.  This  method  unfortunately  produces  a  blurred  image.  To  correct 
this  problem,  a  deblurring  filter  must  be  applied  to  the  data.  It  is  preferred  in  some  applications 
to  filter  the  projections  prior  to  superimposing  (8),  but  the  filtering  can  be  done  on  the  spectral 
data  prior  to  transforming  it.  Since  both  the  frequency  and  spatial  domain  data  are  discrete  in 
nature,  the  location  of  each  individual  image  pixel  may  lie  between  available  data  points.  This 
brings  up  the  requirement  for  interpolation  to  produce  a  more  accurate  image.  The  filtered 
back-projection  method  has  been  shown  to  produce  very  accurate  results  while  also  being 
reasonable  computationally  (1).  Table  1,  taken  from  (1),  compares  computational  speed  of 
the  focusing  methods  described  above. 

1.3  Approach/ Methodology 

Lewitt  clearly  explains  the  development  of  two-dimensional  direct  focus  and  filtered 
back-projection  reconstruction  algorithms  in  (7).  In  this  research  a  three-dimensional  recon¬ 
struction  algorithm  is  developed  using  the  same  type  of  mathematical  development.  Due 
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to  its  computational  efficiency  This  algorithm  is  a  three-dimensional  extension  of  the  two- 
dimensional  Filtered  Back-Projection  algorithm  due  to  its  computational  efficiency  .  This 
algorithm  is  then  applied  to  compact  range  data  taken  on  a  |  scale  model  of  a  C-29  aircraft, 
measured  by  Wright  Laboratory.  A  suitable  method  for  displaying  the  reconstructed  three- 
dimensional  data  is  then  developed  so  that  it  can  be  visually  compared  to  the  actual  physical 
properties  of  the  aircraft. 

1 .4  Summary  of  Chapters 

Chapter  2  covers  imaging  theory  from  its  fundamentals  through  a  detailed  development 
of  three-dimensional  imaging.  Finally  a  discretized  algorithm  is  developed.  The  actual 
application  of  the  theory  is  discussed  in  Chapter  3.  The  target  used  is  described  as  well  as  the 
measured  data  sets.  The  methodology  of  producing  the  actual  images  found  in  this  research  is 
also  discussed.  Three  dimensional  images  as  well  as  two-dimensional  slices  are  included  in 
Chapter  4.  These  results  are  explained  and  compared  to  the  actual  target.  Chapter  5  contains 
conclusions  drawn  from  the  research  and  recommendations  for  future  work. 
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11.  Radar  Imaging  Theory 


2.1  Introduction 

In  this  chapter  radar  imaging  is  presented,  starting  with  its  basic  principles  through 
two-dimensional  Filtered  Back-Projection.  This  two-dimensional  mathematical  basis  is  then 
extended  to  a  three-dimensional  Filtered  Back-Projection  algorithm.  Finally,  this  three- 
dimensional  algorithm  is  then  discretized  for  use  on  a  computer. 


2.2  Two  Dimensional  Radar  Imaging 

Radar  sensors  respond  to  electro-magnetic  waves  which  are  scattered  when  the  propa¬ 
gation  of  incident  waves  is  disturbed  by  the  presence  of  an  object  (8).  When  an  object  such  as 
an  aircraft  is  interrogated  by  a  radar,  certain  features  such  as  wing  tips  and  engines  produce 
“hot  spots”  or  scattering  centers  where  most  of  the  return  energy  is  radiated  from.  A  spatial 
map  of  these  scattering  centers  would  reveal  the  physical  properties  of  the  object.  This  spatial 
distribution  of  reflectivity  corresponding  to  an  object  is  known  as  a  radar  image  (8). 

In  order  to  produce  a  good  quality  radar  image,  a  high  degree  of  both  downrange  and 
crossrange  or  azimuthal  resolution  must  be  achieved.  Figure  1  illustrates  the  definitions  of 
downrange  and  crossrange.  Downrange  location  is  determined  by  measuring  the  time  delay 
between  time  markers  in  the  transmitted  waveform  and  the  received  waveform.  In  pulse  radar 
systems,  these  time  markers  are  the  leading  and  trailing  edges  of  the  pulse.  The  round  trip 
delay ,  r,  for  a  target  at  range,  R,  where  c  is  the  propagation  velocity  is 


So  in  order  to  distinguish  between  two  points  separated  by  a  distance  AR,  you  will 
need  a  pulse  width  of  less  than  T  such  that 


2AR 
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^  Downrange  direction 


Target 

Crossrange  direction 


Radar 

Figure  1.  Definition  of  Crossrange  and  Downrange. 


This  Equation  can  also  be  written  as 


(3) 


or 


(4) 


where,  B  =  is  the  system  bandwidth  as  defined  in  Figure  2.  It  can  be  seen  that  as  the 
system  bandwidth  increases,  the  resolution  becomes  finer,  with  the  best  resolution  occurring 
at  infinite  bandwidth  or  when  the  pulse  is  an  impulse  in  time.  So  in  order  to  get  a  high  degree 
of  downrange  resolution  directly  it  would  be  necessary  to  transmit  the  entire  bandwidth  of 
frequencies  simultaneously  such  as  in  an  Impulse  Radar.  This  would  prove  to  be  extremely 
difficult  if  not  impossible.  Luckily,  the  same  results  can  be  achieved  synthetically  by  making 
several  narrowband  measurements  at  discrete  frequency  increments,  as  shown  in  Figure  2. 
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Af 


Af  2Af  MAf 


Figure  2.  Sampled  Frequency  Spectrum,  B  =  MAf,  1-D  imaging. 

Combining  all  these  narrowband  responses  will  provide  the  same  range  resolution  as  if  the 
entire  frequency  band  was  transmitted  simultaneously  (8). 

Next,  considering  the  crossrange  dimension,  the  radar  aperture  must  be  increased  in 
size  in  order  to  improve  azimuth  resolution.  Ax,  as  shown  in  Equation  5, 


c 


(5) 


where  0  is  the  total  radian  angular  aperture  and  /c  =  /o  +  y^/  as  defined  in  Figure  2. 


Physically  increasing  aperture  size  can  lead  to  the  need  for  extremely  large  antennas. 
A  more  favorable  method  for  our  purposes  is  to  develop  a  synthetic  aperture  by  sequentially 
transmitting  and  receiving  at  discrete  points  along  a  circular  path  to  vary  the  viewing  angle. 
The  received  signals  are  coherently  summed  together  to  produce  signals  equivalent  to  those 
that  would  be  received  by  the  physical  aperture  as  shown  in  Figure  3  (8). 
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Figure  3.  Synthetic  Aperture  Formation,  spotlight  mode  SAR. 

Synthetic  Aperture  Radar  (SAR)  utilizes  both  a  synthetically  expanded  bandwidth  as 
well  as  the  synthetic  aperture  process  explained  above.  SAR  is  utilized  in  many  ways  but 
this  research  is  focused  on  spotlight  Inverse  Synthetic  Aperture  Radar  (ISAR),  in  which  the 
radar  is  stationary  with  the  anterma  fixed  on  the  target,  while  the  viewing  angle  is  changed  by 
rotating  the  target  on  a  pedestal.  It  is  also  assumed  that  the  radar  sensor  uniformly  irradiates 
the  entire  object  as  shown  in  Figure  4. 

As  stated  earlier,  a  radar  image  can  be  looked  at  as  the  sum  of  the  target’s  scattering 
centers  which  can  be  represented  by  the  function  f{x,y).  Our  goal  is  to  determine  the 
locations  of  these  scattering  centers  from  our  reflected  radar  signals  by  producing  an  image. 

First  we  will  define  our  geometry.  As  can  be  seen  in  Figure  4,  the  x  and  y  axis  are 
fixed  with  respect  to  the  target  while  the  u  and  v  axis  are  the  cross  range  and  down  range 
respectively  with  respect  to  the  radar.  For  a  fixed  0,  the  received  signal  for  a  certain  down 
range  point  v  will  be  the  integral  of  the  reflectivity  density  over  u,  assuming  that  all  reflection 
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Angle  of  Rotation 


Radar 

Figure  4.  ISAR  geometry. 

coefficients  have  equal  phase  components.  Note  that  p{v,  0)  in  Equation  6  can  be  looked 
upon  as  a  projection  of  the  reflectivity  on  the  v  axis  at  the  look  angle  0. 
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Figures.  Radar  projection. 

The  total  received  signal  at  a  certain  look  angle  and  frequency,  F(/,  ©),  will  be  the  line 
integral  of  p{v]  0)  along  v  with  a  correcting  round  trip  phase  factor.  Thus 

F{f,Q)=  (  p{v,Q)e~^dv  =  f  p{v,Q)e~^~^dv.  (8) 

oo  */— oo 

Substituting  Equation  6  into  Equation  8  we  obtain 

/°°  r°°  iiTvf 

/  fe{u,v)e  c  dudv  (9) 

-oo  J —oo 


Since  from  Equation  7,  -y  =  a:  sin  0  +  y  cos  0,  F(/,  0)  can  be  expressed  as 


F(/,0)=  r  r 

J  —  OO  J  — OO 

now  defining 

fx  =  f  sin  0 


(10) 


(11) 
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and 


fy^fcosQ  (12) 

F{f,  0)  can  be  rewritten  as 

Fr{fxjy)=f  [  f{x,y)e~^^^^^^^''y^dxdy  (13) 

J— 00  J—oo 

where  Fr{fx,  fy)  is  F{f,  0)  expressed  in  rectangular  coordinates.  Equation  13  is  recognized 
as  the  two-dimensional  Fourier  transform.  Due  to  the  properties  of  the  2-D  Fourier  transform 
we  can  write  (4) 


2-D 

f{x,y)  ^  FrifxJy)  (14) 

2-D 

where  denotes  a  2-D  Fourier  transform  pair.  The  inverse  2-D  Fourier  transform  of 
Equation  13  is  given  by 


f(x,y)  =  r  r  (15) 

J  —00  J—oo 

Thus,  this  equation  establishes  how  a  2-D  image  is  formed. 

In  ISAR,  as  we  defined  it  in  Figure  4,  the  radar  data  is  not  collected  on  a  uniform 
rectangular  grid  since  the  rotation  of  the  target  introduces  polar  aspects.  This  property  would 
result  in  a  distorted  or  “unfocused”  image  if  the  standard  Inverse  Fourier  Transform  is  used, 
especially  when  data  is  collected  over  large  look  angles.  In  order  to  eliminate  this  distortion 
we  must  convert  the  frequency  domain  into  polar  coordinates  as  shown  in  Figure  6.  This 
conversion  is  implemented  by  substituting  Equations  11  and  12  into  Equation  15  as  well  as 
letting 

dfxdfy  =  fdfdQ  (16) 

we  obtain 

f{x,y)=  r  r  (17) 

Jo  Jo 
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Figure  6.  Spatial  frequency  coordinates. 

This  can  also  be  written  as 

f{x,y)=  r  + 

Jo  Jo 

r  r  F(f,  0  +  (18) 

Jo  Jo 

and  since  for  isotropic  scatterers 

F(/,e  +  7r)  =  F(-/,e) 

cos(0  +  tt)  =  —  cos(0)  (19) 

sin(0  +  tt)  =  —  sin(0) 

Equation  17  can  be  rewritten  as 

f{x,y)  =  r  r  F(/,0)e'^(""“®+J'“®®)l/|d/d0.  (20) 

Jo  J—oo 
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Equation  20  produces  a  “focused”  reconstruction  of  the  target.  However  since  the 
standard  rectangular  coordinate  system  is  not  used  the  data  is  not  compatible  with  a  2-D  Fast 
Fourier  Transform  (FFT)  algorithm.  This  makes  computing  a  “focused”  image  as  shown  in 
Equation  20  much  more  computationally  expensive  than  producing  an  “unfocused”  image  (1). 

In  order  to  increase  computational  speed  it  is  also  possible  to  construct  a  “focused” 
image  by  summing,  or  superimposing  one  dimensional  projections  from  all  available  angles 
(8).  This  allows  the  use  of  the  one  dimensional  Inverse  Fourier  Transform  as  follows.  It 
can  be  seen  in  Equation  6  that  the  total  received  signal  at  a  certain  look  angle  is  also  the 
one-dimensional  Fourier  Transform  of  the  projection  p(t;;  ©)  at  that  look  angle.  So,  rewriting 
Equation  20  as 


/•TT  r  /»oo 

f{x,v)=  /  F(/,6), 

JO  iJ —oo 


.^^^(xsin&+y  cos0) 


\f\df 


de 


(21) 


it  can  be  seen  that  the  bracketed  expression  is  a  one  dimensional  Fourier  transform  of  the  total 
received  signal  at  a  fixed  look  angle  0,  multiplied  by  a  filtering  term  |/|.  Let’s  now  call  this 
the  filtered  projection  p{v,  0), 


p{v,e)  =  r  F{f,e)e‘^\fW  (22) 

J  —00 

where  v  =  xsinQ  +  y  cos  0.  Substituting  Equation  22  into  Equation  21  we  obtain 

f{x,y)=  f  p{v,e)de.  (23) 

Jo 

Equation  23  is  known  as  the  Filtered  Back-Projection  method  of  imaging  and  can  utilize  a  one 
dimensional  IFFT  algorithm  to  compute  the  projections.  This  greatly  increases  computational 
speed  (1)  while  also  produces  a  “focused”  image. 
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2.3  Three  Dimensional  Radar  Imaging 

Up  to  this  point  all  discussion  has  been  confined  to  two  dimensions  in  an  effort  to 
develop  a  mathematical  basis  for  three  dimensional  imaging.  In  order  to  obtain  elevation 
resolution  the  target  must  also  be  rotated  in  a  plane  orthogonal  to  the  azimuthal  rotation  as 
shown  in  Figure  7. 


w 


Again  the  u,  v,  and  w  axes  are  fixed  in  relation  to  the  radar  while  the  x,  y,  and  z  axes 
are  fixed  in  relation  to  the  target  such  that  when  ©  =  <^  =  0  the  a:  axis  is  equivalent  to  the  u 
axis,  the  y  axis  is  equivalent  to  the  v  axis,  and  the  z  axis  is  equivalent  to  the  w  axis.  Figure  8 
shows  the  radar  configuration  at  the  0  =  0  =  0  look  angle. 

Assuming  that  the  Fourier  transform  relation  holds  in  three  dimensions 

3-D 

f{x,  y,  z)  ©  fy,  f,)  (24) 


we  can  say 
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i  w 


Figure 8.  ISAR configuration aiO  =  (f)  =  0. 


S(x,v,z)  =  /“  /”  /”  (25) 

y_00  */“00  ^—00 

Again  as  can  be  seen  in  Figure  7,  the  spherical  coordinate  system  is  more  suitable  than 
rectangular  coordinates  for  this  geometry.  Converting  Equation  25  to  spherical  coordinates 
by  the  following  relations  obtained  from  Figures  9  and  10 


fy 

fz 

dfxdfydf. 


f  sin  ©  cos  (j) 
f  cos  0  cos  4> 
fsm(f> 

P  cos  (j)dfdQd(j) 


yields 


(26) 
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Figure  9.  Spatial  frequency  coordinates  in  three  dimensions. 


f{x,y,z)=  r  r 

J —n  Jo  Jo 


A  similar  process  used  to  transform  Equation  17  into  20  will  transform  Equation  27  into 

f{x,y,z)=  r  r  r 
J  — TT  Jo  J  —  00 

(28) 

Again,  to  obtain  the  Filtered  Back-Projection  Equation  we  rewrite  Equation  28  as 

fix,y,z)  =  P  1^°°  ded(f>. 


In  this  case  the  bracketed  expression  is  also  a  one  dimensional  Fourier  transform  at  a  fixed 
azimuth  angle  0  and  elevation  angle  0,  multiplied  by  a  filtering  term,  p  cos  (j).  This  will  be 
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Figure  10.  df^d/ydf^  =  p  cos  (f>dfdQd<j). 


called  a  filtered  projection 

p{v,Q,(f))=  [  F{f,  0,  (j))e^  f  cos  ^df  (30) 

J  — OO 

where  v  =  xsinQ  cos  0  +  y  cos  0  cos  0  +  z  sin  (f).  Substituting  Equation  30  into  Equation  29 
yields 

f{x,y,z)=  f  [  p{v,e,(j))ded(j).  (31) 

J  —It  Jo 

Equation  31  is  known  as  Filtered  Back-Projection  in  three  dimensions.  It  is  this 
equation  that  will  be  utilized  in  the  formation  of  three  dimensional  images  in  this  research. 
But  in  practical  applications  is  is  impossible  to  obtain  an  infinite  bandwidth.  Therefore,  a 
band  limited  version  of  Equation  30  with  frequencies  ranging  from  /o  through  /m  is 


(32) 


pbi{v,  0,0)=  [  F{f,  0,  0)e'^/2  cos  4>df. 

■'fo 

It  is  also  uncommon  to  have  a  complete  set  of  projections  from  every  angle  so  limiting 
Equation  3 1  to  azimuth  angles  from  0o  through  On  and  elevation  angles  from  0o  through  cpi 
gives 


f<t>L  r&N 

fbi{x,y,z)=  /  Pu{v,e,(l))d&d(f).  (33) 

J4>o  JQo 

2.4  Discretized  Three-Dimensional  Algorithm 

In  order  to  implement  this  algorithm  on  a  computer,  a  discrete  version  must  be  developed. 
The  continuous  time  version  already  derived  can  be  sampled  at  discrete  angles  0n,  0  <n  <  N 
and  4)m,  0  <  m  <  M  as  well  as  using  a  sampled  frequency  spectrum  as  shown  in  Figure  2. 
But  inherently,  when  a  signal  is  sampled  in  the  frequency  domain,  aliasing  occurs  in  the 
time  domain  (9).  In  the  downrange  direction  A/  must  be  chosen  such  that  the  downrange 
unambiguous  range,  Rux, 


Fux  — 


C 


(34) 


is  at  least  as  large  as  the  downrange  extent  of  the  target.  The  angular  steps,  A0  and  A0  must 
also  be  chosen  such  that  the  unambiguous  range  in  the  crossrange,  aitd  vertical  direction. 


Ru 


Zi 


Ruy  — 


2/cA0 


(35) 


Ruz  = 


C 

2/cA0 


are  also  greater  than  the  size  of  the  target  in  those  directions. 


(36) 
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Using  the  trapezoidal  rule  (5)  to  discretize  the  filtered  projection  from  Equation  32  we 

obtain 

L  jinvyf, 

Pdi^ki  (f^m)  ~  Fjfli  Qrp  fl  COS  0^.  (37) 

1=0 

Equation  37  can  be  recognized  as  an  IFFT  which  can  be  easily  implemented  using  available 
software. 

In  a  discretized  version  of  Equation  33,  values  must  be  computed  for  each  for  each 
Xq,yri  Zs  pixel  such  that  0  <  q  <  Q,  0  <  r  <  R,  and  0  <  s  <  S.  Since  Vk  is  also  discrete 
it  is  unlikely  that  the  required  equation 


Vk  =  Xq  sin  ©„  cos  (j)m  +  Vr  COS  ©„  COS  (j)m  +  -^r  sin  (f)^  (38) 

will  hold  exactly,  but  rather  the  value  will  fall  between  two  values  such  as  Vk  and  Vk+i. 
Therefore,  in  order  to  obtain  a  value  for  pd  a  linear  interpolation  scheme  must  be  used.  The 
method  used  in  this  research  is  as  follows.  Letting 

Va  =  Xq  sin  ©„  cos  Pm  +  Vr  COS  ©„  COS  pm  +  sin  pm  (39) 


then. 


PiiXat  ©nj  Pm) 


- - —^Pd{Vk,  ©n,  (t>m)  + 

Vk+1  -  Vk\ 


I'U  V  I 


(40) 


so  each  Xq,yr,  pixel  is  computed  as 


AT  M 


fblip^qiyriZ^s)  —  ^Q^p  ^  ^  \Pi{'^ai  ©n;  0m)  • 


(41) 


n=0 1=0 


2.5  Conclusion 

In  this  chapter  radar  imaging  was  presented  from  its  basic  principles  through  two- 
dimensional  Filtered  Back-Projection.  This  two-dimensional  mathematical  basis  was  then  ex¬ 
tended  to  a  three-dimensional  Filtered  Back-Projection  algorithm.  Then,  the  three-dimensional 
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algorithm  was  discretized  for  use  on  a  computer.  The  algorithm  derived  in  this  chapter  is  ap¬ 
plied  to  actual  radar  data  in  the  next  chapter. 
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III.  Application  and  Methodology 


3.1  Introduction 

In  this  chapter  the  target  and  radar  data  used  to  test  the  algorithm  of  Chapter  2  is 
presented.  Then,  the  methods  used  to  generate  results  are  explained  as  are  the  obstacles 
encountered  during  the  implementation  of  these  methods. 

3.2  The  Target 

In  order  to  test  the  algorithm  developed  in  Chapter  2,  three-dimensional  radar  measure¬ 
ments  were  made  on  an  |  scale  model  of  a  C-29  aircraft  in  a  Wright  Laboratory  compact 
range.  This  model  was  accurate  in  every  detail  when  compared  to  the  actual  aircraft  except 
they  put  a  reflective  metallic  film  on  the  windscreen  in  order  to  block  reflections  from  the 
cockpit.  Photographs  and  scale  drawings  of  the  model  are  shown  in  Figures  1 1  and  12. 

3.3  The  Data 

Two  series  of  radar  measurements  were  taken  on  this  target.  One  series  was  taken 
with  azimuth  angles  equaling  —5°  <  ©  <  5°  with  step  size  A0  =  0.04°,  elevation 
angles  equaling  3°  <  ^  <  7°  with  step  size  =  0.04°,  and  frequencies  equaling 
26  Ghz  <  /  <  36  Ghz  with  step  size  A/  =  10  Mhz.  The  transmit  and  receive  polarizations 
were  both  vertical  for  these  measurements.  This  data  set  will  be  referred  to  as  “  the  nose  data”. 
Figure  13  illustrates  this  geometry. 

The  second  series  used  the  same  elevation  angles  and  frequencies  as  the  first  series,  but 
the  azimuth  angles  equaled  —19°  <  ©  <  29°  with  step  size  A©  =  0.04°.  The  transmit 
and  receive  polarizations  were  both  horizontal  for  these  measurements.  This  data  set  will  be 
referred  to  as  “  the  wing  data”.  Figure  14  illustrates  this  geometry. 

Both  data  sets  have  the  same  unambiguous  ranges  as  follows; 
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6.93  meters 


Rux  — 

c 

3x10® 

2/cA0 

-  2(31xl09)(0.04x^) 

— 

c 

2Af 

3x10® 

~  2(10x10®) 

Ruz  = 

c 

3x10® 

2fcA4> 

-  2(31xl09)(0.04Xj|5) 

and  theoretically  provide  resolutions  of, 


15  meters 
6.93  meters 


= 

C 

2/cA0 

ARy  = 

c 

2B 

AR,  = 

c 

2fcA4> 

3x10^ 

2(31xl0«)(10.0Xjfo) 

3xl0» 

2(10x10*) 

3xl0» 

2(31xl09)(4.0x^) 


0.0277  meters 
0.015  meters 
0.0693  meters. 


(42) 


(43) 


3.4  The  Approach 

The  measured  data  or  “raw  data”  is  the  Fp(//,  Qn,4'm)  in  Equation  37  and  as  expected 
the  data  files  are  very  large  at  203.6  megabytes  each.  The  first  step  in  making  an  image  is 
computing  the  filtered  projections,  pd-  This  was  accomplished  using  the  Matlab  program¬ 
ming  language  developed  by  Math  Works  Inc.,  due  to  its  ability  to  quickly  perform  matrix 
computations. 

Unfortunately,  in  spite  of  Equation  43,  it  was  found  to  be  necessary  to  up  sample  the 
data  when  computing  the  IFFT  in  order  to  get  sufficient  resolution  in  the  resulting  image. 
This  is  a  result  of  accumulated  error  from  summing  the  numerous  projections.  This  caused 
a  large  increase  in  data  size  which  required  a  trade-off  between  image  resolution  and  usable 
disk  space.  The  process  of  up  sampling  essentially  provides  additional  interpolated  projection 
values  between  the  originals,  providing  more  accurate  approximations  when  summing  the 
projections  to  compute  a  pixel  value.  It  was  determined,  by  trial  and  error,  that  a  six  times  up 
sample  for  the  nose  data  and  a  seven  times  up  sample  for  the  wing  data  proved  sufficient  to 
provide  an  accurate  image.  In  order  to  demonstrate  the  importance  of  up  sampling  Figure  15 
shows  no  up  sampling,  6  times  up  sampling,  and  20  times  up  sampling  on  one  elevation 
two-dimensional  images  using  the  nose  data.  It  is  seen  that  image  quality  decreases  radically 
as  the  up  sampling  factor  drops  with  the  case  of  no  up  sampling  resulting  in  a  useless  image. 
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However,  the  limited  hard  drive  space  available  for  this  research  resulted  in  the  need  for 
additional  data  reduction.  Since  the  target  is  only  5.2  meters  in  length  and  the  unambiguous 
range,  Ruy,  is  15  meters,  a  50%  data  reduction  is  allowed  by  only  requiring  the  use  of  every 
other  frequency  sample.  Similarly,  since  the  target  is  1.79  meters  long  in  the  ^  direction  a 
A(/)  =  0.08°  still  allows  an  unambiguous  range  greater  than  the  extent  of  the  target  in  that 
direction  as  shown  in  Equation  44. 


Raz  = 


3  X  10® 


2(31  X  109)(0.08  X  jfo) 


=  3.47  meters 


(44) 


This  fact  permits  an  additional  50%  data  reduction  by  making  it  necessary  to  only  use  the  data 
from  every  other  elevation  angle.  The  50  filtered  projection  files,  one  for  each  of  the  elevation 
angles  between  </)  =  3°  and  0  =  7°  would  take  up  approximately  605  megabytes  of  storage 
space  for  the  nose  data  and  705  megabytes  for  the  wing  data.  These  figures  are  still  beyond 
that  which  was  available  for  the  generation  of  this  research,  so  only  37  of  the  elevation  angles 
between  0  =  3°  and  0  =  5.88°  were  used.  This  resulted  in  an  acceptable  new  ^  direction 
resolution  of. 


AR,  = 


3  X  10® 

2(31  X  109)(2.88  X  yfg) 


=  0.0963  meters. 


(45) 


Once  the  filtered  projection  files  were  produced,  this  data  was  transformed,  using  Fortran 
77  code,  to  one  large  three-dimensional  array  according  to  /,  0,  and  0  and  restored  as  an 
unformatted  binary  223 .3  megabyte  file  for  the  nose  data  and  a  260.6  megabyte  file  for  the  wing 
data.  The  use  of  one  large  file  enables  a  great  increase  in  computational  speed  by  allowing 
all  the  data  to  be  loaded  into  resident  memory  under  one  variable,  so  the  calculations  of 
summing  and  interpolation  in  equations  40  and  41  can  be  done  efficiently.  If  enough  resident 
memory  is  available  to  hold  the  filtered  projection  data  as  well  as  the  additional  variables 
to  perform  the  calculations,  an  individual  x,y,z  pixel  can  be  computed  in  0.15  seconds  on  a 
Sun  Microsystems  SPARC  2  workstation.  This  number  greatly  suffers  if  insufficient  resident 
memory  is  available. 
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After  all  the  pixels  are  computed  a  method  of  displaying  the  image  is  needed.  Two 
methods  were  used  in  this  research.  Two-dimensional  slices  were  made  from  the  three- 
dimensional  data  set  using  Matlab  to  show  the  locations  of  the  major  scatterers  and  compare 
them  to  the  actual  target.  Also,  IDL,  developed  by  Research  Systems,  Inc.,  has  a  useful 
volume  visualization  tool  called  SLICER.  This  application  was  used  to  produce  iso-surfaces 
or  three-dimensional  contour  plots  of  a  volume  at  a  given  density  level.  For  example,  in 
medical  imaging  applications,  a  series  of  two-dimensional  images  can  be  created  by  computed 
tomography  or  magnetic  resonance  imaging.  When  stacked,  these  images  create  a  grid  of 
density  measurements  that  can  be  contoured  to  display  the  surface  of  anatomical  structures  (10). 
IDL  also  was  used  to  visualize  slices  of  the  the  volume  data  set. 

S.5  Conclusion 

In  this  chapter  the  |  scale  model  of  the  C-29  and  the  two  series  of  radar  measurements 
used  to  test  the  algorithm  of  Chapter  2  were  presented.  Then,  the  methods  used  to  generate 
results  were  explained  as  was  the  obstacles  encountered  during  the  implementation  of  those 
methods. 
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Figure  11.  Photographs  of  C-29  |  Scale  Model. 
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5.2  meters 


Figure  12.  |  Scale  Dimensions  of  a  C-29. 
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1 .79  meters  5  2  meters 


Azimuth  Angle 

VV  Polarization 
-5°  <  0  <  5° 

A0  =  0.04° 

Elevation  Angle 

3°  <(j)<  7° 

A(/)  =  0.04° 

Frequencies 

26  X  10®  hertz  <  /  <  36  x  10®  hertz 

A/  =  10  X  10®  hertz 

Figure  13.  Geometry  for  Nose  Data  Series. 
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Figure  15.  Illustration  of  effect  of  up  sampling  on  image  quality,  (a)  no  up  sampling;  (b)  6 
times  up  sampling;  (c)  20  times  up  sampling. 
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IV.  Results 


4.1  Introduction 

This  chapter  presents  the  three-dimensional  images  of  the  C-29  aircraft.  These  images 
are  in  the  form  of  iso-surface  images  as  well  as  two-dimensional  slices.  The  images  are 
discussed  and  compared  to  the  physical  properties  of  the  target. 

4.2  Nose  Data 

Three-dimensional  pixel  values  were  computed  using  the  nose  data  set  described  in 
Figure  13  and  the  algorithm  developed  in  Section  2.4.  The  volume  consists  of  720,000  pixels 
and  has  the  following  dimensions, 

—3.0  meters  <  x  <  3.0  meters  stepsize  Ax  =  0.05  meters 

—3.0  meters  <  t/  <  3.0  meters  stepsize  Ay  =  0.05  meters 

—0.5  meters  <  ^  <  2.0  meters  stepsize  Az  =  0.05  meters. 

Figure  16  shows  the  top  view  of  the  iso-surface  for  this  volume.  The  silhouette  of  the  actual 

target  is  overlayed  on  the  image.  The  returns  of  the  engines  and  landing  lights  on  the  leading 
edge  of  the  wings  are  readily  visible.  Figure  17  is  the  same  iso-surface  observed  from  the 
front.  The  vertical  separation  of  the  engines  and  landing  lights  is  evident  and  shows  the 
true  three-dimensionality  of  these  results.  The  pedestal  exhibits  large  returns  due  to  the  fact 
that  VV  polarization  was  used  and  the  pedestal  was  not  subtracted  out  of  the  measurements. 
Sidelobes  from  the  engines  are  also  seen  from  this  view. 

By  looking  at  two-dimensional  slices  taken  at  the  major  scatterers  it  is  easy  to  see  that 
they  line  up  with  the  physical  characteristics  of  the  target.  Figure  18  shows  a  contour  plot 
in  the  xz  plane  at  the  leading  edge  of  the  engines.  An  image  of  this  same  slice  is  shown  in 
Figure  19  from  a  different  aspect.  Figure  20  is  a  contour  plot  in  the  xz  plane  showing  the 
reflections  from  the  APU  intake  and  Figure  21  is  a  slice  taken  of  the  same  area.  Figures  22 
and  23  shows  reflections  from  the  landing  lights. 
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Landing  Lights 
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Cross-Range  (meters) 


Figure  16.  Top  view  of  iso-surface  from  nose 


-3-2-1  O  1  2  3 


Cross-Range  (meters) 

Figure  17.  Front  view  of  iso-surface  from  nose  data. 

Two-dimensional  slices  were  also  made  in  the  xy  plane  to  verify  the  precision  of  the 
return’s  locations  in  the  2:  dimension  when  compared  to  the  target.  Figure  24  is  a  slice  through 
the  engines  in  the  xy  plane  and  clearly  shows  the  engine’s  reflections.  Figure  25  is  a  slice 
through  the  location  of  the  APU  intake  in  the  xy  plane  and  Figure  26  shows  a  slice  through 
the  landing  lights. 

Figures  27  through  33  are  slices  taken  in  the  xz  plane  of  the  whole  volume  computed 
from  the  nose  data.  These  slices  are  taken  every  0.1513  meters  going  from  front  to  back. 
Figures  34  through  36  are  slices  taken  in  the  xy  plane  of  the  same  volume.  These  slices  are 
taken  every  0.1530  meters  from  top  to  bottom.  Figures  37  through  42  are  expanded  views  of 
some  of  the  slices  taken  at  critical  areas.  They  are  also  presented  with  the  outline  of  the  plane 
overlayed  on  them. 
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Vertical  [z  axis]  (meters)  Vertical  [z  axis]  (meters) 


Cross-Range  [x  axis]  (meters) 


Figure  18.  Contour  plot  in  xz  plane  of  engine  returns. 
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I 

I  Slice  at  y=-0.0756 

L _ _ _ _ 
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Down-Range  [y  axis]  (meters) 

Figure  19.  Slice  in  xz  plane  of  engine  returns. 
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Figure  20.  Contour  plot  in  xz  plane  of  intake  returns. 
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Figure  22.  Contour  plot  in  xz  plane  of  landing  lights  returns. 
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Down-Range  [y  axis]  (meters) 

Figure  23.  Slice  in  xz  plane  of  landing  lights  returns. 
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Figure  24.  Slice  in  xy  plane  of  engine  returns. 
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Figure  25.  Slice  in  xy  plane  of  APU  intake  returns. 
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Cross-Range  [x  axis]  (meters) 


Figure  26.  Slice  in  xy  plane  of  landing  lights  returns. 


VerticaJ  [z  axis]  (meters)  Vertical  [z  axis]  (meters) 


Figure  27.  Nose  Data:  xz  plane  slices  at  (a)  y  =  3.0000  meters  (b)  ?/  =  2.8487  meters  (c) 
y  =  2.6975  meters  (d)  y  =  2.5462  meters  (e)  y  =  2.3950  meters  (f)y  =  2.2437 
meters. 
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meters. 
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Figure  29.  Nose  Data:  xz  plane  slices  at  (a)  y  =  1.1849  meters  (h)y  =  1.0336  meters  (c) 
y  =  0.8824  meters  (d)  y  =  0.7311  meters  (e)  y  =  0.5798  meters  (f)y  =  0.4286 
meters. 
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Vertical  (zaxbl  (meters)  Verticat  (z  axis)  (meters)  Vertical  (z  axis)  (meters) 
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(c)  (d) 


(e)  (f) 

Figure  30.  Nose  Data:  xz  plane  slices  at  (a)  y  =  0.2773  meters  Qo)  y  =  0.1261  meters 
(c)  y  =  —0.0252  meters  (d)  y  =  —0.1765  meters  (e)  y  =  —0.3277  meters  (f) 
y  —  —0.4790  meters. 
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Cross-Range  [x  axis]  (meters) 


Cross-Range  [x  axis]  (meters) 


(e)  (f) 

Figure  31.  Nose  Data:  xz  plane  slices  at  (a)  y  —  —0.6303  meters  (b)  y  =  —0.7815  meters 
(c)  y  =  —0.9328  meters  (d)  y  —  —1.0840  meters  (e)  y  =  —1.2353  meters  (f) 
y  =  —1.3866  meters. 
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Cross-Range  [x  axis]  (meters) 


Cross-Range  [x  axis]  (meters) 


Figure  32.  Nose  Data:  xz  plane  slices  at  (a)  y  =  —1.5378  meters  (h)y  =  —1.6891  meters 
(c)  y  =  —1.8403  meters  (d)  y  —  —1.9916  meters  (e)  y  =  —2.1429  meters  (f) 
y  =  —2.2941  meters. 
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Vertica/  [z  azisl  (metere)  Vertical  fz  axis]  (melers) 


Cross-Range  (x  axis]  (meters) 


Cross-Range  [x  axis]  (meters) 


Nose  Data:  xz  plane  slices  at  (a)  y  =  —2.4454  meters  (b)  y 
(c)  y  =  —2.7479  meters  (d)  y  =  —2.8992  meters. 


2.5966  meters 


Down-Range  ^  axis]  (meters)  Down-Range  ^  axis]  (n>6(ers)  Down-Range  (y  axis]  (meters) 


Figure  34.  Nose  Data:  xy  plane  slices  at  (a)  z  —  2.0000  meters  (b)  2  =  1.8469  meters  (c) 
2;  =  1.6939  meters  (d)  z  =  1.5408  meters  (e)  ^  =  1.3878  meters  (f)  -2  =  1.2347 
meters. 
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Figure  35.  Nose  Data:  xy  plane  slices  at  (a)  z  =  1.0816  meters  (b)  z  =  0.9286  meters  (c) 
2  =  0.7755  meters  (d)  z  =  0.6224  meters  (e)  z  =  0.4694  meters  (f)  z  =  0.3163 


meters. 
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Figure  36.  Nose  Data:  xy  plane  slices  at  (a)  z  =  0.1633  meters  O^)  -2^  =  0.0102  meters  (c) 
z  =  —0.1429  meters  (d)  2  =  —0.2959  meters  (e)  ^  =  —0.4490  meters. 
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Figure  37.  Nose  Data:  Slice  in  xz  plane  of  landing  lights  scattering. 
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Figure  38.  Nose  Data:  Slice  in  xz  plane  of  APU  intake  scattering, 


Vertical  [2  axis]  (meters) 
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Figure  39.  Nose  Data:  Slice  in  xz  plane  of  engines  scattering. 
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Figure  40.  Nose  Data:  Slice  in  xy  plane  of  landing  lights  scattering. 
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Figure  41.  Nose  Data;Slice  in  xy  plane  of  engines  scattering. 


57 


xy  slice  at  z=1 .2347 
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Figure  42.  Nose  Data:  Slice  in  xy  plane  of  APU  intake  scattering 


43  Wing  Data 

The  wing  data  described  in  Figure  14  provides  the  most  dramatic  results.  The  HH 
polarization  makes  the  pedestal  less  of  a  reflector.  The  19°  to  29°  azimuth  angles  are  normal 
to  the  leading  edge  of  the  right  wing  and  leading  edge  of  the  right  horizontal  stabilizer.  This 
angular  range  allows  the  wing  and  the  horizontal  stabilizer  to  be  the  major  scatterers.  Also, 
this  makes  it  easy  to  see  the  vertical  separation  between  the  wing  and  the  tail. 

A  volume  consisting  of  864,000  pixels  was  computed  over  the  following  dimensions, 

—3.0  meters  <  x  <  3.0  meters  stepsize  Ax  =  0.05  meters 
—3.0  meters  <  y  <  3.0  meters  stepsize  Ay  =  0.05  meters 
—0.5  meters  <  z  <  2.5  meters  stepsize  Az  =  0.05  meters. 

Figure  43  is  a  top  view  of  the  iso-surface  for  this  volume  and  shows  reflections  from 
the  wing  and  horizontal  stabilizer.  Figure  44  shows  the  vertical  separation  of  the  wing  and  the 
horizontal  stabilizer  coincides  with  their  physical  location  on  the  target.  The  reflections  at  the 
bottom  of  Figure  44  are  caused  by  the  pedestal. 

Figures  45  through  5 1  are  slices  taken  in  the  xz  plane  of  the  entire  volume  computed 
from  the  wing  data.  These  slices  are  taken  every  0.1513  meters  going  from  front  to  back. 
Figures  52  through  55  are  slices  taken  in  the  xy  plane  of  the  same  volume.  These  slices  are 
taken  every  0. 1525  meters  from  top  to  bottom.  Figures  56  and  57  are  expanded  views  of  slices 
in  the  xy  plane  of  the  horizontal  stabilizer  and  wing  respectivly. 

4.4  Conclusion 

Images  computed  by  utilizing  the  three-dimensional  algorithm  of  Chapter  2  were  pre¬ 
sented  in  this  chapter.  The  accuracy  of  this  algorithm  was  verified  by  comparing  iso-surface 
images  as  well  as  two-dimensional  slices  to  the  physical  properties  of  the  target.  The  images 
demonstrated  that  scattering  primarily  occurs  at  physical  discontinuities  on  the  target.  In 
the  case  of  the  nose  data,  VV  polarization  was  used  causing  surfaces  with  vertical  edges, 
such  as  the  pedistal,  to  be  dominant  scatterers.  Since  HH  polarization  was  used  for  the  wing 
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Cross-Range  [x  axis]  (meters) 

Figure  43.  Top  view  of  iso-surface  from  wing  data. 
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Figure  44.  Front  view  of  iso-surface  from  wing  data. 

data,  surfaces  with  horizontal  edges  are  the  dominant  scatterers.  This  combined  with  azimuth 
angles  normal  to  the  wing  and  horizontal  stabilizer’s  leading  edge  cause  these  “edges”  to  be 
the  dominant  scatterers. 
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Figure  45.  Wing  Data:  xz  plane  slices  at  (a)  y  =  3.0000  meters  (b)y  =  2.8487  meters  (c) 


y  =  2.6975  meters  (d)  y  =  2.5462  meters  (e)  y  =  2.3950  meters  (f)  y  =  2.2437 
meters. 
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Figure  46.  Wing  Data:  xz  plane  slices  at  (a)  y  =  2.0924  meters  (b)  y  —  1.9412  meters  (c) 
y  =  1.7899  meters  (d)  y  =  1.6387  meters  (e)  y  =  1.4874  meters  {f)y  =  1.3361 
meters. 
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(f) 


Figure  47.  Wing  Data:  xz  plane  slices  at  (a)  y  =  1.1849  meters  (b)y  =  1.0336  meters  (c) 
y  =  0.8824  meters  (d)  y  =  0.7311  meters  (e)  y  =  0.5798  meters  (f)y  =  0.4286 
meters. 
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Vertical  [z  axis]  (meters)  Vertical  [z  axis]  (meters)  Vertical  [z  axis]  (meters) 


Figure  48.  Wing  Data:  xz  plane  slices  at  (a)  y  =  0.2773  meters  (b)  t/  =  0.1261  meters 
(c)  y  =  —0.0252  meters  (d)  y  =  —0.1765  meters  (e)  y  =  —0.3277  meters  (f) 
y  =  —0.4790  meters. 
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Vertic^  (z  axisj  (meters)  Vertical  (z  axis]  (meters)  Vertical  (z  axis)  (meters) 


Cross-Range  [x  axis]  (meters) 


Cross-Range  [x  axis]  (meters) 


-3-2-10123-3-2-10123 
Cross-Range  (x  axis]  (meters)  Cross-Range  [x  axis]  (meters) 


(e) 


(f) 


Figure  49.  Wing  Data:  xz  plane  slices  at  (a)  y  =  —0.6303  meters  (b)  y  =  —0.7815  meters 
(c)  y  =  —0.9328  meters  (d)  y  =  —1.0840  meters  (e)  y  =  —1.2353  meters  (f) 
y  =  —1.3866  meters. 
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Vertical  [z  axis]  (meters)  Vertical  [z  axis)  (meters)  Vertical  (z  axis]  (meters) 


-3-2-10123-3-2-10123 


Cross-Range  (x  axis]  (meters) 


Cross-Range  [x  axis]  (meters) 


-1  0  1 
Cross-Range  (x  axis]  (meters) 


-1  0  1 
Cross-Range  [x  axis]  (meters) 


-3-2-10123-3-2-10123 
Cross-Range  [x  axis]  (meters)  Cross-Range  [x  axis]  (meters) 

(e)  (f) 

Figure  50.  Wing  Data:  xz  plane  slices  at  (a)  y  —  —1.5378  meters  (b)  t/  =  —1.6891  meters 
(c)  y  =  —1.8403  meters  (d)  y  =  —1.9916  meters  (e)  y  =  —2.1429  meters  (f) 
y  =  —2.2941  meters. 
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Vertical  (z  axis]  (meters)  Vertical  [z  axis]  (meters) 


Figure  51.  Wing  Data:  xz  plane  slices  at  (a)  y  =  —2.4454  meters  (b)  y  =  —2.5966  meters 
(c)  y  =  —2.7479  meters  (d)  y  =  —2.8992  meters. 
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Figure  52. 


Wing  Data:  xy  plane  slices  at  (a)  z  =  2.5000  meters  (b)  2:  =  2.3475  meters  (c) 
2:  =  2.1949  meters  (d)  z  =  2.0424  meters  (e)  2:  =  1.8898  meters  (f)  ^  =  1.7373 


Down-Range  [y  axis]  (meters)  DoMwRange  ly  axb]  (meters)  Down-Range  [y  axis]  (meters) 


-2-1  0  1  2  3 

Cross-Range  [x  axis]  (meters) 

(b) 


-2-10123 
Cross-Range  [x  axis]  (meters) 
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Cross-Range  [x  axis)  (meters) 


Cross-Range  [x  axis]  (meters) 

(c) 


Cross-Range  [x  axis)  (meters) 


Figure  53.  Wing  Data:  xy  plane  slices  at  (a)  =  1.5847  meters  (b)  2:  =  1.4322  meters  (c) 

2;  =  1.2797  meters  (d)  2  =  1.1271  meters  (e)  2  =  0.9746  meters  (f)  2  =  0.8220 
meters. 
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Range  ly  axis)  (meters)  Down-Range  |y  axis)  (meters) 


-3-2-1  0  1  2  3 

Cross-Range  (x  axis]  (meters) 


-3-2-10123 
Cross-Range  [x  axis]  (meters) 


Figure  54.  Wing  Data:  xy  plane  slices  at  (a)  2  =  0.6695  meters  (b)  z  =  0.5169  meters 
(c)  2  =  0.3644  meters  (d)  2  =  0.2119  meters  (e)  2  =  0.0593  meters  (f) 
2  =  —0.0932  meters. 
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-Range  [y  axis]  (meters) 


-2-10123 
Cross-Range  [x  axis)  (meters) 
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-2-10123 
Cross-Range  [x  axis]  (meters) 
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Figure  55.  Wing  Data:  xy  plane  slices  at  (a)  z  =  —0.2458  meters  (b)  2  =  —0.3983  meters. 

xy  slice  at  z=2.1949  meters 


-2-1012 
Cross-Range  [x  axis]  (meters) 


Figure  56.  Wing  Data;  Slice  in  xy  plane  of  horizontal  stabilizer  scattering. 
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Down-Range  [y  axis]  (meters) 


xy  slice  at  z=1 .2797  meters 


3 

-3  -2-1  0  1 

Cross-Range  [x  axis]  (meters) 


Figure  57.  Wing  Data:  Slice  in  xy  plane  of  wing  scattering. 
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V.  Conclusions 


5.1  Summary  of  Work 

In  Chapter  2,  basic  radar  imaging  theory  was  presented  from  its  basic  principles  through 
two-dimensional  Filtered  Back-Projection.  With  this  basis  in  hand,  three-dimensional  Filtered 
Back-Projection  was  then  developed  and  translated  into  a  discrete  algorithm.  Chapter  3 
presented  the  C-29  one-third  scale  model  used  to  generate  data  for  this  research  as  well 
as  two  series  of  three-dimensional  radar  measurements  taken  by  Wright  Laboratory.  The 
methodology  used  to  convert  these  radar  measurements  into  actual  radar  images  was  recounted 
while  some  of  the  obstacles  confronted  in  this  process  were  touched  on.  The  results  were 
presented  in  Chapter  4.  Iso-surface  images  were  presented  from  different  view  angles  and 
two-dimensional  slices  were  included  of  the  major  scattering  features. 

5.2  Discussion  of  Results 

This  research  demonstrates  the  ability  to  produce  accurate  three-dimensional  images 
given  three-dimensional  radar  data.  The  major  problem  confronted  in  this  research  was  the 
lack  of  memory  available  in  terms  of  both  hard  drive  space  and  RAM.  This  problem  limited  the 
amount  of  up  sampling  allowed  when  computing  the  projections  and  ultimately  led  to  inferior 
resolution  as  well  as  spurious  responses  in  the  final  image.  Theoretically,  sufficient  resolution 
should  have  been  available  without  upsampling  the  original  data  as  shown  in  Equation  43. 
However,  due  to  accumulated  error  caused  by  summing  the  numerous  projections,  large 
errors  resulted  in  the  final  image.  In  order  to  reduce  this  effect,  up  sampling  as  well  as  the 
interpolation  step  described  in  Equation  40  was  necessary  (1).  It  was  determined  through  trail 
and  error  that  at  the  least,  up  sampling  by  a  factor  of  6  for  the  nose  data  and  7  for  the  wing  data 
produced  a  minimumly  acceptable  image.  Using  a  factor  below  those  given  above  resulted 
in  poor  image  quality  as  shown  in  Figure  15.  Ideally,  up  sampling  by  a  factor  of  20,  as  was 
used  in  Anderson’s  two-dimensional  research  (1),  would  produce  excellent  image  quality  but 
would  result  in  a  dataset  far  to  large  to  work  with  for  three-dimensional  imaging.  However, 


74 


it  is  still  obvious  from  the  images  shown  in  Chapter  4  that  this  three-dimensional  algorithm 
produces  accurate  scatterer  location  in  the  2:  dimension.  If  sufficient  memory  were  available 
to  up  sample  at  a  higher  rate  much  better  resolution  can  be  achieved  in  the  resulting  image. 
The  most  dramatic  results  were  produced  using  the  wing  data.  It  is  clear  by  viewing  Figure  44 
that  a  true  three  dimensional  representation  of  the  target  was  computed. 

5.3  Recommendations  for  Future  Research 

Primarily,  any  future  work  in  the  area  of  three-dimensional  imaging  should  be  focused 
on  reducing  the  memory  requirements  of  computing  an  image.  Inherently  the  “raw”  three- 
dimensional  radar  data  will  require  extremely  large  storage  space.  The  aforementioned  need 
to  up  sample  the  projection  data  causes  a  multitudinous  increase  in  storage  space  requirements 
over  this  original  data.  It  would  be  of  great  benefit  to  look  into  an  alternative  method 
of  interpolation  other  than  up  sampling.  One  possibility  is  illustrated  in  Oppenheim  and 
Shafer  (9).  The  process  of  up  sampling  is  equivalent  to  applying  an  ideal  low-pass  filter  in 
the  frequency  domain  or  convolving  with  a  sine  function  in  the  time  domain.  The  use  of  a 
convolution  to  determine  the  interpolated  values  eliminates  the  need  to  up  sample.  But  care 
should  be  taken  not  hinder  the  computational  speed  of  computing  a  pixel  in  the  process.  The 
fact  that  such  a  high  number  of  pixels  need  to  be  computed  to  obtain  a  useful  image,  for  example 
864,000  pixels  for  the  nose  data  in  this  research,  even  a  small  decrease  in  computational  speed 
could  cause  a  large  increase  in  the  time  required  to  compute  a  volume  of  pixels. 

Additional  research  can  possibly  be  done  in  the  area  of  artificially  extending  the  band¬ 
width  of  the  radar  data  through  the  use  of  parametric  methods.  This  process  would  effectively 
remove  any  noise  in  the  computed  projections  and  possibly  produce  a  clearer  image  in  the 
end. 
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Appendix  A.  Computer  Code 


This  appendix  contains  the  computer  code  used  to  compute  the  images  in  this  research. 
Section  A.l  contains  the  code  that  computes  the  filtered  projections  and  stores  them  according 
to  elevation  angle  in  Matlab  .mat  files.  Section  A.3  converts  the  .mat  files  to  an  3-D  unformatted 
fortran  array.  Section  A.4  contains  the  code  that  actually  computes  the  3-D  volume  from  the 
filtered  projections. 


A.l  projection.m 

% 

%  File:  projection. m 

%  Author:  Jack  D.  Pullis 

%  Purpose:  Compute  filtered  projections  on  raw  radar  data  and  store 

%  according  to  elevation  angle. 

\ 

% 

% 

f=  [26002300:20000:36002300]*10*3; 

ph=[4. 92: .08:6.04] ; 

o-3e8; 
m=7 ; 

df=abs(f(2)-f(l)); 
dth=abs(th(2) -th(l) ) ; 
dph'»abs(ph(2)  -ph(l) ) ; 

for  el=l ; length (ph) ; 

E=num2str( (30-ph(el) )*100 ) ; 
file=['C29_3D_hhN.O' ,E] ; 
fid=fopen(file) ; 

buf=fread(fid, [2008,251] , 'float32'); 
fclose( 'all' ) 

G=buf (7:2:2008, : )+j *buf ( 8 : 2 : 2008,  : ); 

[M,V]=size(G); 

M=501; 

Ru=0.5*o/df;  %  calculate  the  unambiguous  range 

dr=Ru/(m*M);  *  calculate  the  range  resolution 

r=(0:m*M-l)*dr;  »  form  the  downrange  vector 


filt=cos(ph(el)*(180/pi) )*( (f ( : ) .*2)*ones(l,V) ) ;  %  calculate  imaging  filter 

corr=exp( j*4*pi*f (l)*r'/c)*ones(l,V) ;  %  shift  to  baseband  for  IFFT 

rc=corr . *ifft{G(l : 2 : 1001, : ) . *filt,m*M) ;  *  compute  filtered  projections 

for  ci=l:V; 

rc{ : ,ci)=fftshift(rc( : ,ci) ) ; 
end 

E2=num2str(ph(el) ) ; 

file2=[ '/horae/maroonil/c29/C_29_3D_hhNel' ,E2] 

eval(['save  '  file2  '  rc'])  %  store  projections  in  elevation  file 


%  speed  of  light 
%  oversample  factor 
%  calculate  delta-f 
%  calculate  delta-theta 
\  calculate  delta-phi 
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end 


A.2 

0 
c 
c 
c 
c 
c 
c 
c 

c--- 
character*50 

coraplex*8  rc(3006,251,37) 
coraplex*16  r(754506) 
real*4  ph(34) 

integer*4  fp,matOpen,p 
integer  ran,  stat 
integer*4  raatGetMatrix 
integer  mxGetM,  raxGetn,  raatClose 

c . . . . . . 

fill='/liorae/raaroonil/c29/Binproj/datahn2  .all' 

open ( unit=55 , file=fill , f orm= ' unformated ' , status- ' unknown ' ) 

ph(l)-3.00 
do  20  1-2,37 
ph(I)=ph(I-l)+.08 

20  continue  ' 


do  10  1-1,37 

open (unit-52,  status- ' scratch ' ) 

write(52,ll)  '/home/raarconil/c29/C_29_3D_hhNel' ,  ph(I) 

rewind (unit-52) 
read (52, 12)  filel 

11  format(A32,F4.2) 

12  format (A36) 
close(52) 

write (6,*)  filel 

fp  -  matOpen (filel,  'r') 

p  -  raatGetMatrix (fp,  'rc') 

ran  =  mxGetM(p)  *  mxGetN(p) 

call  mxCopyPtrToCoraplexl6(mxGetPr(p) ,  mxGetPi(p),  r,  ran) 
stat  =  matClose(fp) 
call  mxFreeMatrix(p) 

write (6,*)  ran 
do  30  K=l,251 
do  40  J-1,3006 

indx-indx+1 

rc(J,K,I)-r(indx) 


convert./ 

file:  convert. f 
author:  Jack  D.  Pullis 

purpose:  Convert  .mat  files  containing  filtered  projections 
to  a  unformated  fortran  3-D  array  then  save  it 
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40  continue 
30  continue 

indx=0 


10  continue 

write{55)  ro 
close(55) 

end 


A.3  imageSd.f 

c 

c  file:  iraage3d.f 

C  author:  Jack  D.  Pullis 

0  purpose:  Compute  3-D  Filtered  Back -Projection  volume  from 

c  pre-stored  projection  file  named  data. all.  The  resulting 

0  data  is  stored  in  a  Matlab  .mat  file  under  the  name 

0  rdata.mat 

0 
c 
o 

0 . 

dimension  ph{101) ,th(251) , f (501) ,st(251) 

dimension  x(90) ,y(90) ,z{5) , v(251) ,ct(521) ,sph(101) ,cph(101) 

real  *4  ph 

integer  ran,  stat,inc,inol,rl,inc2,Q,indx,13 
integer  f p , raatOpen , p , mxCreat ePul 1 
complex  *8  ro(3507, 251, 37) ,oum 
complex  *16  data(40500) 
integer  matClose 
characterise  filel,file2 

complex*8  rcl ( 3006 ,251,19), ro2 ( 3006 ,251,18) 
character*50  fill,fil2,fil3 

c . . . — . 

c 

c  Load  Projection  data 

0 

c - - - 

fill= ' /home/fishl/j  saoohin/data . all ' 

open ( unit»55 , file»fill , form- ' unformated ' , status- ' unknown ' ) 

read(55)  rcl 

close(55) 

c . . . 

c 

c  Load  parameter  data 

o 

0 . . . . 

c=3e8 

open ( unit=l , file- ' outtest . dat ' , status- '  old ' ) 
read(l,*)  ph,th,f,x,y,z 

ml»6 

nth-251 
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nph=101 

M=501 

nx=90 

ny=90 

nz=l 

df=abs(f(2)-f(l)) 
dth=abs(th(2)-th(l)) 
dph=abs(ph(2)-ph(l) ) 

Ru=0 . 5*c/df 
dr=Ru/(ral*M) 

do  20  1=1, nth 

radi=th(I)*{3. 14/180) 

st(I)=sin(radi) 

ct(I)=cos(radi) 

20  continue 

do  30  1=1, nph 

sph(I)=sin(ph{I)* (3. 14/180)) 
oph(I)=oos(ph(I)*{3 .14/180) ) 
30  continue 


0 . 

c 

0  Compute  each  pixel 

0 

c . 

do  40  zi=l,nz 
do  50  xi=l,nx 
do  60  yi=l,ny 
oura=(0.0  ,  0.0) 
inc“0 

do  70  1=1,73,2 
inc=inc+l 
zsoal=sph(l) 

do  75  12=1, nth 
xsoal»st(12)*oph(l) 
ysoal=ct ( 12 ) *oph ( 1 ) 


v(  12 )=x(xi) *xscal+y (yi ) *ysoal+z ( zi ) *zscal 
75  continue 


do  76  13=1, nth 

rp=v(13)/dr+l+(M*ml)/2 

rl=int{rp) 

oura=oum+(rl+l-rp)*ro(rl,13,inc)+(rp-rl)*ro(rl+l,13,inc) 
76  continue 


70  continue 

data(yi+( (xi-l)*ny)+{zi-l)*ny*nx)-oum 
60  continue 
write{6,*)  xi,zi 
50  continue 
40  continue 
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fp  =  matOpen( 'rdata.mat' ,  'w' ) 
p  =  inxCreateFull(nx*ny,  nz,  1) 

call  mxCopyComplexl6ToPtr( data,  mxGetPr(p),  nixGetPi(p),  nz*nx*ny) 

call  itixSetName(p,  'R') 

stat  =  matPutMatrix(fp,  p) 

stat  =  matClose(fp) 

call  raxFreeMatrix(p) 


end 
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